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Abstract 

We study the basic features of the two-dimensional quantum Hubbard Model at 
half-filling by means of the Liischer algorithm and the algorithm based on direct 
update of the determinant of the fermionic matrix. We implement the Liischer idea 
employing the transfer matrix formalism which allows to formulate the problem on 
the lattice in (2 + 1) dimensions. We discuss the numerical complexity of the Liischer 
technique, systematic errors introduced by polynomial approximation and introduce 
some improvements which reduce long autocorrelations. In particular we show that 
preconditioning of the fermionic matrix speeds up the algorithm and extends the 
available range of parameters. We investigate the magnetic and the one-particle 
properties of the Hubbard Model at half-filling and show that they are in qualitative 
agreement with the existing Monte Carlo data and the mean-field predictions. 
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9. 



1 Introduction 



It is very well known that quark degrees of freedom are difficult to include in numer- 
ical simulations of quantum systems [Q . Although the substantial progress has been 
made in this field the search for new algorithmic techniques is an active research 
area. At present probably the best method is to integrate the quark degrees of 
freedom and simulate the theory with effective action. However such action usually 
contains determinant which is a nonlocal object and its evaluation is very expen- 
sive in computer time. Recently a new computational scheme has been proposed 
by Liischer who showed, in the QCD case, how to approximate determinant of the 
Dirac operator by a local bosonic action 0. So far there is no complete discussion 
of properties of this algorithm and it is of great interest to perform some tests of its 
performance. 

Besides QCD there is a wide class of physically interesting problems where effec- 
tive technique for simulating fermionic degrees of freedom is of crucial importance. 
The Hubbard Model is an example of a conceptually simple model with effec- 
tive Coulomb repulsion between electrons which for many years has been the basic 
framework for studying strongly correlated electrons in solid state physics. The in- 
terest of the Hubbard Model has been renewed after finding that it can qualitatively 
describe the magnetic properties of the high Tc superconducting materials above 
the temperature of the superconducting transition 0]. There is also a conjecture 
suggested by a variety of calculations that when two dimensional system is doped 
slightly away from half-filling the repulsion between electrons can give rise to su- 
perconductivity Until now this hypothesis has not been confirmed or rejected. 
The Hubbard Model is very well suited to tests of the new fermionic algorithms be- 
cause it contains all details of fermionic formulation and describes very interesting, 
and rather unexplored from numerical point of view, physics of strongly correlated 
electrons. 

In this paper we apply the Liischer idea to the Hubbard Model. We show how 
to study basic features of strongly correlated electrons and note about the problems 
which can be solved by numerical simulations. We study the efficiency of the Liischer 
algorithm in comparison with the algorithm based on direct update of the deter- 
minant of the fermionic matrix. We propose also some improvements which reduce 
long autocorrelation times. Some conclusions from this dicussion can be applied in 
implementation of the Liischer algorithm in other areas. 

The paper is organized as follows. Section 2 is assigned to the review of the 
basic properties of the Hubbard Model. We emphasize some important not fully 
solved problems. Section 3 covers path integral formulation of the Hubbard Model. 
In section 4 we review fermion Quantum Monte Carlo methods and introduce two 
new algorithms. One is the implementation of the Liischer method to the Hubbard 
Model. The second is a modification of the algorithm based on direct update of the 
determinant of the fermionic matrix. In section 5 we compare dynamical properties 
of these two algorithms and present some physical results obtained from simulations 
of half-filled band Hubbard Model. The comparison of these results with previous 



findings and theoretical predictions gives credence to our methods of simulations. 
Section 6 contains a summary and our conclusions. 

2 The Hubbard Model 

The Hubbard Model is defined by the Hubbard Hamiltonian 




where Ojo- is the annihilation operator of an electron with spin a located at site 
i of d dimensional lattice. In this paper however we will concentrate mainly on 
the case d = 2, which is physically the most interesting. The summation in the 
kinetic term runs over nearest neighbors (each symmetric pair of {i,j) is counted 
twice) and the hopping energy K can be taken as a unit of energy. The interaction 
term takes into account short-range effective Coulomb repulsion. The term with 
the chemical potential determines the band filling. The half-filled band has /i = 0. 
There are three dimensionless parameters which determine the behaviour of the 
system, interaction energy U/K, band filling ^/K and temperature T/K. In this 
parameter space the Hubbard Model exhibits various phases. 

Originally the Hubbard Model was introduced as the simplest system which 
might exhibit an insulating (Mott) state. It is exactly soluble only in one dimension. 
In this case the integral equation for the density of states can be solved in a closed 
form . It follows that the system is insulating for any nonzero U and conducting for 
[7 = 0. Thus there is no Mott transition at finite U . However even for Id system we 
have only partial information. For example, the correlation functions are not known. 
In 2 and higher dimensions many different approximate techniques have been used. 
Mean field approximation [0 , variational techniques |^ and other were successful in 
describing the properties of strongly correlated electrons. They are very important 
because they give insight into physics. Unfortunately it is very difficult to judge the 
accuracy of approximate techniques beyond the perturbative limit. Therefore it is 
not surprising that they often give conflicting results. Numerical simulations can fill 
this gap and deliver nonperturbative results for comparison. 

Some properties of the Hubbard Model on square lattice may be understood 
from simple qualitative considerations. Let us start from the noninteracting limit. 
The dispersion relation 

Cfc = — 2_ft'(cos /c^ + cos /Cy), (2) 

establishes the band structure. The total bandwidth Ei, = 8K. The ground state 
is obtained by populating the states from the bottom of the energy scale up to the 
Fermi level. In the half-filled case the rectangular shape of the Fermi surface is 
determined from the condition ei? = as shown in Fig. 1 . 

The existence of parallel sectors on the Fermi surface which are separated by 
the vector Quest = (iTT, tt) (nesting) has important physical consequences and is 
crucial for the developing of the magnetic instabilities. Indeed the weak interaction 
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could induce the process of exchanging electrons with opposite spins from opposite 
sides of the Fermi surface. It involves however the exchange of momentum Quest 
which is characteristic for antiferromagnetism. Obviously without detailed theory 
it is impossible to make final conclusions. For example, there is also a process of 
quasi-forward scattering with the exchange of momentum (0, 0) which would give 
rise to ferromagnetism and from simple qualitative arguments it is impossible to 
judge which is dominant. 

On the other hand, in the strong coupling limit {U —>■ oo) the interaction term 
forces electrons not to reside on the doubly occupied site. The hopping term can 
be treated as a perturbation and the ordinary expansion in the parameter K/U is 
applicable. The effective Hamiltonian is the Quantum Heisenberg Antiferromagnet 

'^0 = —j^^SiSj, (3) 

^ {id) 

with the exchange coupling 2K'^/U. The above result is valid for the half- filled 
system in any dimension and lattice. 

Therefore it is plausible that the Hubbard Model may possess ground state with 
nontrivial magnetic properties. More quantitative description can be obtained by the 
mean field theory. One assumes that the ground state has some sort of magnetic 
order and the problem can be solved by means of variational principle. Detailed 
calculations for the antiferromagnetic ground state is given in the Appendix A. 
Here we stress main results. One possible solution is the ferromagnetic state which 
can exists for sufficiently large U . The critical value is given by Stoner's criterion 

Ucritpiep) > 1 (4) 

with piep) = J2k^{^F — Cfe) being the density of states on the Fermi surface. But 
there is also an antiferromagnetic solution for any U with the smaller energy which 
is a better candidate for the true ground state of the half-filled Hubbard Model. Nu- 
merical simulations fully confirm that expectation and the antiferromagnetic order 
has been explicitly proven for [/ < 2 |p. It is also believed that antiferromag- 
netism exists for any nonzero U. Mean field solution gives also new band structure 
Ek = ±-^e| -|- with the gap 2 A above the Fermi surface characteristic for an 
insulator. 

A lot of single-particle properties can be obtained from the study of the one 
particle Green's function defined by 

G,.,(r,t;r',t') = -z(Ta,(r, t)a^,(r', t')) (5) 

where acr{r,t) are annihilation operators in Heisenberg picture. T orders product of 
operators at earlier time to the right of operators at later time, keeping track of the 
sign changes associated with Fermi statistic; i. e. for two fermionic operators A{t) 
and B{t) we have 

TA{t)B{t') = A{t)B{t')e{t - t') - B{t')A{t)e{t' - t). (6) 



For the noninteracting system, which is translationally invariant and time inde- 
pendent the Fourier transform of G{r, t; r', t') is a function of 2 variables 



The poles of the propagator exhibit typical one particle spectrum. In the weakly 
interacting case we can expect that physical low-energy excitations look like weakly 
interacting fermions. The propagator retains its pole structure with a renormalized 
dispersion relation eren{k). Thus it is expected that the propagator near the Fermi 
surface should look like 

lim ^(fc.Cj) = lim ( ^- + Greg] (8) 

with a certain regular part Greg- The wave function renormalization Z measures the 
discontinuity of the occupation number Uka = iGaaik.t = 0~) at the Fermi surface. 
The system which exhibits such behavior is called a Fermi liquid. 

The ground state of the half-filled Hubbard Model has long range antiferro- 
magnetic order being the result of strong interactions. Thus we may expect quite 
different behaviour than the mentioned above. In fact the renormalization constant 
Z vanish and simulations fully confirm the lack of sharp Fermi surface [|lT| |jl2[ . 



What happens after doping still remains unclear. The approximate analytic 
methods are not conclusive and one has to rely on numerical calculations. Unfortu- 
nately in this case the fermionic determinants are not positive definite. As will be 
discussed in Section 4 this difficulty called the sign problem significantly restricts 
the range of predictions. Nevertheless quantum MC can still deliver useful informa- 
tions. For example, it was confirmed that antiferromagnetic order is destroyed even 
for small doping [|TT| At the same time however it is not certain weather the 
ground state becomes the Fermi liquid. The simulations performed at the quarter 
filling strongly suggest the existence of the sharp Fermi surface. The situation is 
unclear at small doping where the sign problem is very severe. Recent simulations 
[0 suggest the vanishing Z as the size of lattice increases at U = 4. However in 
the last calculations the small magnetic field was introduced to alleviate the sign 
problem and it is not clear how it can infiuence results. Thus more work should be 
done to clarify this issues. Should the non Fermi liquid behaviour be confirmed the 
holes would not propagate as free entities. One possible scenario is binding of holes 
and possible existence of superconducting phase. 

The problem of the existence of the superconducting state has been studied 
directly in numerical simulations [|rT| []T^. The appropriate observable is the corre- 



lation function for pairs 

P(r) = (A1a,+,). (9) 
The annihilation operator for pair is defined as 

Aj = ai|(ai+£; + tti-xi ± ai+yi ± ai+yi) (10) 



where the (+) sign correspond to extended s* wave, and the (— ) sign to d^z-y^ wave, 
as can be seen from rotational symmetry |]14|. In the superconducting state these 



correlations should diverge with the spatial lattice volume. A widely used quantity 
to monitor such behaviour is the susceptibility 

X = T.Pir). (11) 

r 

Current simulations indicate only very weak dependence of the pair correlations on 
lattice size what does not confirm the existence of the superconducting phase. One 
may argue however that the lattice sizes available for computation are too small 
while comparing with the correlation length for the pair. 

The Hubbard Model can be modified in many ways. Adding the next-nearest 
neighbour hopping has very important consequences because the Fermi surface 
ceases to be nested. In this case critical value of U exists for the appearance of 
ferromagnetism. There is also some evidence that such model would be more rel- 
evant to the problem of high Tc superconductivity and recent numerical results 
indicate enhancement of correlation functions for pairs These issues are still 
intensively studied. 

Contrary to the results discussed above the superconducting phase has been 
found for the attractive Hubbard Model i.e. with the reversed sign in front of U 



[[Tql . In this case the QMC simulations are much simpler since the sign problem 
does not occur. It might appear that physics of the attractive and the repulsive 
model should be very similar because both models are related by a change of sign 
in U and the exchange of chemical potential and magnetic field. However there is 
no evidence that the attractive model can serve us as an effective model for some 
range of parameters of the repulsive model, where one would expect pairing. Thus, 
at present simulations of the attractive model simply are not a solution of the sign 
problem. 



3 Path integral formulation of the Hubbard Model 

In this section we derive the path integral formulation of the Hubbard Model. De- 
coupling of interaction terms by means of the Hubbard-Stratonovich transformation 
and applying transfer matrix formalism allows to write the expression for the parti- 
tion function in such a form that it can be used as a starting point for Liischer local 
bosonic approximation. 

Our derivation of the Euclidean path integral partition function follows closely 



that presented by Creutz All arguments apply essentially for both: the at- 

tractive and the repulsive Hubbard Model. However there are also some important 
differences which will be stressed. We begin with the Hamiltonian Eq. (|l]), which 
may be rewritten as 

T-(- = -K ^aflja - - ^^lY - H ^i'T- (12) 

{i,j)a i i(T 
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We assume that i and j are sites on a two dimensional square lattice with = 
sites. This formula differs from the standard one for the repulsive Hubbard Model 
by an irrelevant additive constant. 

All thermodynamic functions for a many particle system in the temperature T 
can be obtained from the partition function 



Z = Tr[exp-{n/T)]. 



(13) 



Tr denotes the sum over the complete set of physical states. The thermal expectation 
value of some physical observable, say O, is defined by 



(O) =Tr[Oexp-(n/T)]/Z. 



(14) 



Unfortunately, it is usually very difficult to evaluate the trace in Eq. (0) exactly for 
physically interesting systems. Lattice field theory offers an interesting possibility 
to evaluate the expectation values for finite systems employing stochastic methods. 

In the lattice formulation one interprets the inverse temperature, 1/T, (which is 
also referred to as Euclidean time) as the extra dimension and discretizes it dividing 
the interval (0,/5 = 1/T) into Nt slices of length e = j3/Nt. The partition function 
can be written as 

Z = Tr[{e-^^/''^f\ (15) 

In the next step the separation between kinetic /C and interaction V term is per- 
formed according to the Trotter formula 



(16) 



which starts the series of approximations in which the terms of order of 0{{(3 /N-i)'^) 
are consistently neglected 



-m/Nt 



exp 



exp 



K(3 

UI3 
2N, 



E 



— all a 



(17) 



As all the terms presented in the last exponent commute with each other we can 
treat them as ordinary numbers. In particular one can use the identity 

^-axV2+bx ^ (27^)^/2e''/2^ a > , (18) 

which is a continuous version of the Hubbard- Stratonovich transformation. With 
this in mind we can introduce the integration over auxiliary set of variables Ai 
located on the lattice sites to decouple the quartic term 



-m/Nt 



exp 



(27r) 



(19) 
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Note an ambiguity which appears in the last expression, namely which sign 
should have square root of (riji — rijj^. Two possible options are equally good and 
we arbitrary choose one. In principle the asymmetry between up and down com- 
ponents should vanish after performing functional integral over Ai and should not 
affect the measured observables. However due to further finite Nt approximations 
the corresponding results for two spin components will differ in real simulations. 

To clarify the positivity of the final fermionic determinant it is convenient to 
formally recover the symmetry between up and down polarizations by performing 
the particle-hole transformation. Introducing the creation operators for hole 6| = 
(^—ly^+'^y and denoting simply by Oj we get 



exp 



L * (id) 

j:{U(3/Nty/'A{ala, + bjh - 1) - ^ E(«l«^ + ^Ih) 



(20) 



which is explicitly symmetric in Oj and bi operators. The sum in the exponents can 
be expressed as products with help of the following relations 



exp {alttj) = 1 -|- alttj {i ^ j) , exp (aajai) = 1 + a\ai{e" — 1) 



(21) 



After collecting time slices together and shifting the integration the resulting 
trace of normal ordered products takes the form: 

Z = exp [7V2/3(c//2) + /,)](27r)-^'^*/' 

!\DA,^,]e-^^,^^y^TTY[, : [n(,,)[l + {K(3 /N,)a\a,] 

n.(l + ala,{exp [{U P / N^f^ A,^t - {U - /.)/5/iV,] - 1}) 

({a,at,/.}^ {6, 6t, -;,})] : (22) 

For every time slice, labeled by a discrete index t = 1, . . . ,Nt the anticommuta- 
tion relations can be used to convert this expression into a sum of normal ordered 
operators, neglecting higher orders in (3/Nt . The trace can be expressed as a Grass- 
mann integral. In Appendix B we present the general formalism of the Grassmann 
variables and the derivation of the partition function. The final result is 



Z^H j [£>A,t]e" det M+ det M_, 



(23) 



where M is some normalization constant. The matrices M_|_ and M_ are specified 
by their elements between arbitrary vectors ^0* and 



Kf3 



+ Y.^lti^i,t{^^V [{up/Ntf'^A^t - ([/ T /x)/?/Ar,] - 1}. (24) 



The antiperiodic boundary conditions ipi^ = —ipi,Nt taken in the t direction 
and reflect the Grassmann nature of the fermionic fields (see Appendix B). The M 
matrices are not Hermitian. It seems to be an intrinsic property of the systems with 
Galilean group of symmetry and can not be simply avoided ||18[ . This fact is rather 
disadvantageous for effective implementation of the Liischer method as can be seen 
later. 

Our expression for the partition function is rather nonstandard. We decided 
to work in the representation with dimension V = NN^ to avoid multiplication of 
matrices present in the standard formulation (cf Eq. (^Q])). Moreover to obtain the 
simple structure for the matrix M we performed the series of approximations in 
the parameter e = (3/Nt treating the higher order terms in such a way to get the 
simplest possible form of matrices M±. Due to these approximation our expression 
for the partition function becomes exact only as Nt ^ oo. In principle all measured 
quantities should be obtained by extrapolation to this limit. 

The half-filled case is favoured from the computational point of view because 
then the statistical weights are positive 

det M+ det M_ = det (25) 

since M = M_ = M+. The same it is not true outside half- filling and the sign 
problem occurs (see Section 4). Attractive Hubbard model is more tractable as the 
distribution is positive definite for any fi. The interaction term can be written as 

-^EK+^a)'- (26) 

and the derivation follows essentially the same steps. However in contrast to the pre- 
vious case after performing the Hubbard-Stratonovich transformation the transfer 
matrix remains symmetric in the spin components and the particle hole transforma- 
tion is superfluous. The statistical weights in the final result are evidently positive 
definite 

Z = W j[DAi^t]e~^^'^'^'^'''\deiM+f. (27) 
The matrix M_|_ is still given by Eq. (p^) with U replaced by \U\. 

4 Algorithms 
4.1 Introduction 

The numerical difficulties with fermionic fields come from their anticommuting na- 
ture. The transfer matrix is an operator in a Grassmann space and its matrix 
elements cannot be interpreted directly as a probability in Monte Carlo algorithms. 
For many interesting problems the formalism similar to the described above can 
be used to integrate the fermionic degrees of freedom. Once the final expression 



in 



for the partition function is given in closed numerical form standard techniques of 
stochastic sampling can be apphed. 

Unfortunately fermionic determinants appearing in the partition function are 
non-local objects. This means that local updates require calculations which depend 
on the state of the whole system. Thus to update the fields individually on each 
lattice site, or link, soon becomes prohibitively expensive in computer time as the 
size of system grows. At present, for example in the QCD case, the best known 
algorithms update many degrees of freedom simultaneously. 

Additionally to the nonlocality of the action one has often another difficulty 
with lattice fermions known as a doubler problem. It is for example evident in the 
lattice formulation of QCD where the species of quarks are doubled in the simplest 
scheme incorporating the chiral symmetry. There are two popular approaches for 
dealing with this difficulty. In the Kogut-Suskind formulation each site carries only 
a single component of the Dirac spinor. In the Wilson formulation chiral symmetry 
of Lagrangian is explicitly broken and it is expected to be restored only in the 
continuum limit |T^. The Hubbard Model is physically defined on the lattice and 
hence the doubler problem does not appear in its path integral formulation. 

In this section we review some fermionic algorithms. We discuss the algorithms 
based on evolution equations which are widely used in QCD and algorithms devel- 
oped specially for simulations of the Hubbard Model. Afterwards, we present two 
particular algorithms which we use for simulations of the Hubbard Model. One is 
the application of the Liischer idea to the path integral representation, which we 
discussed in the previous section. The second is a modification of the algorithm 
based on direct computation of the determinant and was introduced as a reference 
point for the comparison of results and efficiency of the Liischer algorithm. 

4.2 Algorithms based on evolution equations 

To be more specific we concentrate on the action typical for fermionic problems 



We denote the bosonic part of the action depending on the field A (which is the 
Hubbard-Stratonovich field for the Hubbard Model) by Sb{A). The matrix Ai 
contains all details of fermionic formulation and its explicit form depends on the 
model. We assume that matrix Ai is positive definite. For example in simulations 
of the half-filled Hubbard Model M = M with matrix M given by Eq. (H). 

The determinant of the matrix M. can be expressed as a gaussian integral over 
a set of auxiliary complex fields 



The problem of evaluating determinant is now reduced to the of inversion of matrix 
M.. Each time components of the A field are changed the evaluation of new Ai^^ 




(28) 




(29) 
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matrix is required. The problem is slightly simplified because for local actions A4 
is sparse. In this case there are some useful iterative schemes appropriate for this 
task. Commonly used is the conjugate gradient algorithm [pJl . 

However the inversion of the matrix still requires a lot of CPU time. Hence to 
reduce the computational effort many approximate algorithms have been proposed. 
One of the simplest and interesting is the pseudofermion method [^]. In this ap- 
proach one considers only small changes in the A field linearizing the change of the 
action 

^J^^^^.TriM-^'-^) (30) 

dA dA \ dA J 

with 

S = SB-Tr\ogM{A). (31) 

For the local matrix Ai the quantity ^ can be easily calculated. The elements of 
the /iA~^ can be obtained as the appropriate correlation functions 

{M-% = {x)x^). (32) 

The expectation value is taken over complex fields x? called pseudofermions which 
are distributed according to the formula 

P(X) cxexp(-x*-Mx)- (33) 

Fields X can be easily simulated giving the estimation for the elements of matrix 
A4~^. Such estimation is usually done once per full sweep of the variables A. The 
assumption of small changes, which can be realized by taking sufficiently small 
step for proposing trial values of the A field is very important. The algorithm is 
only approximate and in principle all measured quantities should be obtained by 
extrapolation with the size of the step going to zero. In practice such analysis is 
rarely made due to the insufficient computing resources. The algorithm described 
above usually has long autocorrelations and is inefficient. However it was the starting 
point for developing new better approaches. 

Another very friutfull idea refers to stochastic equations, which can be used to 
generate the fields with the probability distribution ([28|). For sake of simplicity we 
restrict the discussion to a single degree of freedom. Further generalization to the 
field theory is straightforward. Let us consider a particle with mass m moving in 
a potential V{x). To the Newton's equation of motion we add also a drag force 
proportional to the velocity and a randomizing force 

d^x dV dx f 2a\^^'^ , , 
m— = -— - a— + — r]iT) (34) 
ar ax ar \ 7 / 

with a and 7 being free parameters. The random force is white noise {'r]{T)r]{T') = 
6{t — t')). Detailed definition of the distribution p{ri) of will be given later. The 
coefficient in front of 77 is a matter of convention. After introducing the momentum 
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p the equation can be rewritten as a system of two first order equations 

dp dV ap /2a\^^'^ . , 

dr dx m \ 'y J 

dx p 
dr m 

For the simulation purposes we divide the time r into shces of size e. The quantities 
x' and p' at time r + e can be calculated from x and p at time r from the formulas 



(35) 



P 



X 



p + e 



dV ap f 2a \ ^^'^ , , 

-+[ — ] V{r) 

dx m V 7 / 



= x + 



ep 



m 



(36) 



To minimize the errors of discretization special form of update was introduced. The 
momentum p is calculated first and this new value is used to update x. In the limit of 

deterministic evolution this "leap-frog" scheme of integration preserves phase space 
volumes ( the Jacobian of the transformation equals 1) and is reversible under the 
change of the sign of e . These properties are very useful in making the algorithm 
exact. 

The quantity r] is generated, independently for each discretized slice of the evo- 
lution time, according to distribution 



Pin) - A/^e 



-ri'e/2 



(37) 



which implies that 



p{rj)rj^drj 



1, 

0, 





1 

2 



(38) 



0(e-^/2) j>3 



Let us denote a probability density of finding the state with given x and p by 
P{x,p). Updating the states gives a new ensemble with the probability distribution 



P'{x',p') 



J dx dp P{x,p)P{x,p ^ x' ,p') 
dx dp drj p{r])P{x,p) 



xS \ p' — p — e 



dV ap 
dx m 



1/2 



2a 
7 



xd(x - X ). 

m 



A little algebra yields the result 
P'{x,p) = P{x,p) 

+ € 



'dHdP dHdP" 



dx dp dp dx 



fld'^P pdP P" 
y7 9p2 m dp m 



(39) 



+ 0(6^), (40) 



whrere H is the Hamiltonian corresponding to the original Newton's equation of 
motion H = + V{x). Eq. pO|) is a Fokker-Planck equation for the evolution 
of the probability distribution P{x,p). It is now readily verified to order e that 
stationary distribution for such evolution is the simple Boltzman weight 

P{x,p) = exp{—jH{p,x)} (41) 

which can be obtained through the very long evolution. 

In the limit m —>■ the Eq. (|34D leads to the Langevin equation in the rescaled 




time (r — ar) 

^ = -^+(-1 Vir). (42) 

In the limit a ^ one recovers the deterministic limit called also microcanonical. 
Since in that case we get back to the Newton's equation and the energy remains 
unchanged during evolution the temperature should be measured by some sort of 
thermometer. In our example it could be the average kinetic energy of the particle 
i/cT = (^)- To change the temperature one should start from different initial 
configuration. 

Both approches: Langevin and microcanonical were used in numerical simula- 
tions of fermions p2| [|23|. In the standard version Hybrid Monte Carlo (HMC) 



algorithm p4[ makes use of deterministic evolution equations. For illustration pur- 
poses we consider the HMC algorithm for the half-filled Hubbard Model. Introducing 
for every Ai the corresponding momentum pi the Hamiltonian of the system can be 
written as 

HiA,p) = ^ Ep' + ^ E 4' + E ^KM^M)-]^,, (43) 

i i i,j 

where the indices i and j run over the whole space-time lattice. Initially the gaussian 
momenta p are chosen and fields are generated from gaussian vectors r by = 
M'^^Ajr. Then the system evolves for Nmic steps. The evolution equations are: 

A . dH 



The integration of Eqs. (H3) is performed according to the leap-frog scheme with a 



discrete step e, Nmic^- is the trajectory length. The conjugate gradient algorithm is 
required to compute the new vector (M'^M)~^0 in each step of the evolution. 

Due to the finite step errors the energy does not remain constant during the evo- 
lution. Thus to ensure the exactness of the algorithm additional global reject/accept 
step is needed. According to Metropolis scheme one can fully restore the detailed 
balance by accepting the new configuration (p', A') with the probability 

Pace = min[l, exp {H{p, A) - H{p\ A'))] (45) 

The acceptance rate behaves as 

Pace = erf c{cN^,ct^^), (46) 
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where V is the total volume of the lattice and c is a constant factor. The errror 
function erfc is defined as 



erfc{x) = I e"*' dt . (47) 

y/n Jx 

If the trajectory length is fixed the acceptance falls exponentially with e^. Thus to 
minimize the autocorrelations one should keep the step size as large a possible still 
maintaining reasonable acceptance. 

Currently it is the most widely used algorithm in simulations with dynamical 
quarks in QCD. In the optimal situation its numerical complexity behaves with the 
volume of the system as VV^^'^, as can be readily verified from Eq. (|46|) which is 
only slightly worse than the linear dependence for bosonic algorithms. However 
these algorithms are still rather complicated comparing with the simplicity of pure 
bosonic algorithms and suffer from the strong autocorrelations between generated 
configurations. The proof of the detailed balance of the HMC algorithm requires 
exact reversibility and it is unclear how finite numerical accuracy of calculations can 
influence the reversibility of the algorithm and further the physical results |2^. As 



can be seen later Liischer method has no such difficulties. The resulting action is 
local and systematic errors are under control. 



4.3 Algorithms for simulation of the Hubbard Model 

It is rather surprising that the algorithms described above have not been widely 
used for simulations of the Hubbard Model. Only recently extensive studies of the 
attractive Hubbard Model has been performed with the help of HMC algorithm |]2B 



Today's leading algorithm [|TT| is based rather on the exact numerical evaluation of 
the determinants appearing in the expression for the partition function. Assuming 
that the two dimensional lattice with spatial volume = N^Ny and Nt time slices 



is considered, the partition function in Hirsch formulation reads |27 



Z = detM+(s)detM-(s). (48) 

Here 

M- = I + B%B%^_,...B'[ (49) 

and 

S± = e^^-'^(')e-^"^ (50) 

I is the unit matrix, = 5i,jSi,i and K is the matrix representation of the 

kinetic part of Hamiltonian. The s is the Ising like Hubbard-Stratonovich field 
which has two values +1 or —1. The index i labels the lattice sites and / the 
time slices. The M" matrices are N x N dimensional. A direct evaluation of the 
determinant requires 0{{NNtY) operations. Fortunately due to the very simple form 
of matrices there is a method of updating the determinant with 0{N^Nt) operations 
needed to perform full sweep through the lattice. The early implementations of this 
algorithm suffered from numerical instabilities especially strong in the limit of low 



temperatures and extrapolation to the ground state was impossible. It turned out 
that numerical instabilities occurred during multiplication of the badly conditioned 
matrices B with almost linearly dependent columns. A simple remedy was the 
periodic reorthogonalization of the fermionic matrices and gave very good results. 
At present simulations of the half-filled Hubbard Model can deliver results with 
decent precision in the wide range of parameters. 

An extremely difficult and not fully solved problem, referred to as the sign prob- 
lem, concerns the simulation of the fermionic system when the determinant is not 
always positive. This is the case of the Hubbard Model outside half-filled band. A 
simple method for dealing with this difficulty consists of attaching the unwanted sign 
of the probability measure to the observable. More precisely, if P is not positively 
definite then the expectation value of some observable O can be written as 

{0)p = -r-- — r5V\ — 

{sign{P))iPi 

where (. . .)\p\ means the expectation value with the respect to the absolute value 
of P. Unfortunately the above equation, exact in principle, dramatically increases 
statistical errors when the expectation value of sign becomes small. It is remarkable 
that the sign problem is especially strong just below half-filling where there is a 
greatest chance of finding the superconducting phase . 



It seems that modification of the standard Quantum Monte Carlo algorithms 
called Projector Quantum Monte Carlo will be more useful for obtaining the ground 
state property of the Hubbard Model. In this approach, proposed in Ref. 



for bosonic systems and applied further to the Hubbard Model the projected 
partition function is introduced 

Q = {^T\exp[-/3nMT) (52) 

\iPt) is a trial wave function nonorthogonal to the ground state (usually the Slater 
determinant). The Hamiltonian is used to systematically project the trial wave 
function onto the ground state. Thus in the limit of large imaginary time (3 the 
properties of the ground state can be obtained. The quantity Q can be evaluated 
by the Monte Carlo methods after rewriting it as a path integral. 

The very interesting point is the possibility of circumventing the sign problem in 
this approach. It can be achieved by the appropriate choice of the trial function. It 
can be also handled by applying another statistical weight which gives the identical 
distribution in the low temperature limit with the reasonable expectation value of 
sign. The construction of better probabilistic weights is rather problem dependent, 
some simple examples are given in Ref. p9[. Recently similar philosophy has been 



applied to the noninteracting electron system with chemical potential [^. The de- 
sired distribution is constructed iteratively with the help of the cluster algorithm. It 
was shown that for this simple model the sign problem can be completely eliminated. 
Of course it would be of great interest to extend this method to more complicated 
models. 

Efficiency of the Projector Quantum Monte Carlo algorithm depends on the 
type of the ground state for noninteracting system (U=0). Usually if the ground 



state is unique the convergence is easily achieved. Conversely in the case of the 
degenerate ground state the algorithm suffers from the negative sign problem and 
poor statistics. 

According to Eq. the allowed states in the momentum space can be divided 
into groups (shells) with a given energy. Degeneracy in the ground state appear as a 
result of the existence of not completely occupied shell. More precisely, if m electrons 
with the specified polarization occupy the shell which can contain maximally up to 
n such electrons it leads to a degeneracy 



N'i={ 1\- (53) 




Although there is no analytical expression which gives the degeneracy as a function 
of lattice size and total number of electrons it is easy to compute this degeneracy 
numerically dividing the states in the momentum space into shells. As an example 
we consider a system with almost half-filled band, 7 electrons up and 7 electrons 
down on the 4x4 lattice. The ground state has a degeneracy of Nd = 29, as far 
as the sub-space of total momentum Q = is considered. Therefore we may expect 
rather poor efficiency of PQMC algorithm in this physically important case. 

Quite different numerical approach is based on the exact diagonalization of many 
particle Hamiltonian. The main problem is the enormous dimension of the Fock 
space Dpock = 4^. Fortunately symmetries of the problem often significantly reduce 
this number. Although the dimension of matrices precludes their direct diagonaliza- 
tion it is relatively easy to obtain eigenvectors for few largest and smallest eigenvalues 



using Lanczos method ||31|. They contain all necessary information needed to com- 
pute physical observables. Additional very nice feature of this method is the ability 
to obtain the interesting quantities directly in the Minkowski time (real frequency 
properties). In Quantum Monte Carlo calculations are carried out in imaginary 
time and the statistical errors inherent to the Monte Carlo make the analytical 
continuation difficult. 



4.4 Local bosonic algorithm 

Recent proposal of Liischer has attracted a lot of interest p|. Originally being 
introduced as an alternative to the Hybrid Monte Carlo algorithms in QCD, soon 
has been apphed in QMC simulations of the Hubbard Model The basic idea 
of Liischer consists in approximating of the inverse of the fermionic matrix with 
a polynomial P„ of even degree. Its particular form being proposed is built from 
Chebyshev polynomials and gives rapid uniform convergence to the function 1/x in 
the interval (e, 1), with e being a small positive number 

lim PJx) = (54) 
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The relative error of the approximation defined as Rn{x) 
tially bounded 



\RJx)\ <2 



1 + 



n+l 



xPn{x) — 1 is exponen- 



(55) 



Approximation region should cover the whole spectrum of the fermionic matrix or 
at least its significant part. Therefore it is impossible to apply directly the original 
Liischer polynomial P„ to the nonhermitian matrix M with eigenvalues distributed 
on the complex plane. Instead, we wrote the polynomial approximation for the 
matrix Q'^Q = M"^ M/ X^ax, where Xmax is the largest eigenvalue of M'^M. Hermitian 
and positive definite matrix Q^Q has properly distributed eigenvalues in the interval 
(0,1). 



1 



2n 



k=l 



(56) 



k=l 



since the roots = + iPk of real polynomial come in complex conjugate 
pairs. Their values are analytically known with the accuracy of 0{l/n). More 
precise values can be simply obtained by applying few Newton-Raphson iterations. 
After introducing auxiliary fields (pk the local bosonic representation for the partition 
function of the Hubbard Model reads 



J [dAd(f)]t 



i i,j k=l 



where 



with M given by Eq. 
lattice with volume V 
of the real and the imaginary part of 



{Q^Q-akf + f3l 



(57) 



(5J 



.Here the indices i and j run over the whole space-time 
N"^ X Nt- The action can be rewritten also directly in terms 



^ = E4V2 + E E E 



4j 



i,j r=0,l k=l 



(59) 



Some remarks are in order. The resulting action 5* contains only local inter- 
actions. Hence we can study the fermionic system using bosonic Monte Carlo 
techniques. However some complications are introduced by the presence of next- 
to-nearest neighbors interactions. The condition number of matrix Q^Q is square 
of condition number for M and we may generally expect the poorer behaviour of 
numerical procedures. Moreover it complicates the implementation of this algorithm 
on vector and parallel computers. The analogical action for QCD is simpler because 
one make use of additional symmetry of the Dirac operator. 
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To define the normalized matrix Q we used the upper bound for the Xmax which 
can be easily derived 

Xma. < 2[A^„,(AU) + max(D)2] (60) 

where D is diagonal part of M = A + D. Large magnitudes of the A fields are 
naturally cut by the gaussian part of the action and we may safely assume maximal 
value for the diagonal elements; e. g. the A[x) < 4 would correspond to four standard 
deviations assuming that the gaussian part of the action dominates. 



4.5 Exact algorithm 

Below we discribe algorithm based on direct computation of fermionic determinant. 
We use this algorithm as a kind of reference point which, first allows us by compar- 
ison to control the accuracy of the Liischer algorithm, and second gives us a chance 
to compare an efficiency. 

The idea of exact algorithm is based on the updating scheme allowing to com- 
pute changes in matrix while making trial changes in matrix M. In order to 
update the field A^ located on the site k one must calculate the ratio of the fermion 
determinants. Since the only A dependent elements of matrix M lie on the diagonal 
we consider the following change in the matrix M 

M' = M + A, (61) 

where A is matrix with one nonzero element, say k-th on the diagonal, Ajj = SikSkjd. 
Then 

= det (/ + A) = 1 + M,-,M (62) 

is determined completely by elements of matrix M~^. Once the trial change has been 
accepted the updated can be evaluated from the Shermann-Morrison formula 

[0 

M'- = M- - ^i:^. (63) 

This process is economical from the computational point of view since one up- 
date of requires OiV"^) operations comparing with OiV^) operations needed 
to evaluate the determinant with the brute force method. Today's leading QMC 
simulations of the Hubbard Model are essentially based on it [|lTl . As has been pre- 
viously reported numerical instabilities often appear in such calculations especially 
at low temperatures. However our particular formulation of path integrals does not 
require multiplication of badly conditioned matrices and we believe that it is free 
of this difficulty. Indeed we performed thousands sweeps at (3 as large as 8 without 
accumulating numerical errors. However the additional cost comes from working 
with larger matrices. 



IP 



5 Results 



In this section we present implementation details of the Liischer algorithm. We 
compare the accuracy and efficiency of the local bosonic algorithm with the exact 
algorithm described in the previous section. We use these algorithms to study basic 
features of the Hubbard Model. In our studies we looked mainly at the magnetic 
properties of the Hubbard Model and one-particle Green's functions (shape of the 
Fermi surface and effective hopping). Due to the sign problem our discussion is 
restricted to the half-filled band. 



5.1 Implementation details and the errors of the polynomial 
approximation 



The goal of the MC simulation is to generate A and (p fields with the disribution 
To achieve this goal we alternately update all components of A and all components 
of 0. We ffist consider the update of fields. Because of the gaussian form of 
dependent part of the action in Eq. (0) we applied the heat bath algorithm 



to update this sector. This is implemented as follows. For every (f)^^] component 
we compute the conditional parameters of its gaussian distribution while keeping 
other components fixed. Then the value obtained from gaussian random number 
generator is assigned to the (j)^^]. In a full sweep through the lattice these operations 
are performed for each i, r and k. 

The effective action depends on A in more complicated way and and the heat- 
bath algorithm is not longer useful. Instead we use a more general Metropolis 



algorithm The trial part of a single update relies on proposing a change 

A^ = Ai + a{x - 0.5). (64) 

The Metropolis step accepts the proposed change with the probability 
min(e~^'^, 1), where A^* is corresponding change of the action. This guarantees 
the detailed balance condition, a; is a random number uniformly distributed in the 
interval (0, 1). The parameter a introduces the scale of changes and it is very impor- 
tant to tune it to get minimal autocorrelation times. Usually it is choosen to keep 
the acceptance ratio near 50% . It is also very important to maintain the balance 
between updates of A and fields. Earlier studies of QCD implementation sug- 
gest that one should perform updates of the A field much more frequently than the 
fields to minimize autocorrelation times [^]. In our simulations ten Metropolis 
updates of A field follow one heat-bath generation of fields. 

Since the CPU cost is expected to be at least proportional to the number of fields 
n it is important to adjust this number carefully. One sees from the Eq. ( pSj) that to 
decrease e one has to increase n to keep the error of approximation constant. In the 
case of the approximation of the determinant of the matrix Q'^Q one can easily guess 
that the most economic choice of e should be comparable to the smallest eigenvalue 
of the matrix Q'^Q. To make this discussion more precise we introduce as a measure 
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\R2n{x)\ 


10-2 


10-^ 


10-4 


e = 0.003 


26 


34 


45 


e = 0.001 


41 


60 


78 


e = 0.0005 


59 


84 


110 



Table 1: The number of fields needed to achieve different levels of relative errors |i?2n(2;)| 
on the interval (e, 1) 



of the error of the polynomial approximation the quantity 

5=|yV^-l| (65) 

with 

y = det{Q^QPniQ^Q)) (66) 

The power 1/V properly normalizes the quantity y per one degree of freedom. Fig. 
2 shows 5 as a function of e for few numbers of fields measured on one typical 
configuration and confirms simple expectations [Xmin = 0.0032 ). With this hint 
we estimated the number of fields n to be about 50-100. This guaranteed that the 
relative error |-R2n| was smaller than 10-^ in the whole interval (e, 1). If one accepts 
bigger relative errors, for example 10-^ these numbers can be reduced but in this 
case we found discrepancies between exact algorithm and Liischer implementation 
(Table. 1). 



5.2 Dynamical properties of the Liischer algorithm and its 
modifications 

In order to establish the computational effort required to generate one uncorre- 
lated (independent) statistical event we studied the integrated autocorrelation times 
for different observables. Let {Ot}(t = 1, . . . ,Ns) be an ordered sequence of data 
produced by our algorithm with the mean value /x. We define the unnormalized 
autocorrelation function 

C(t) = (0,0,+i)-/i2, (67) 
normalized autocorrelation function 



Pit) = C(t)/C(0), (68) 
and integrated autocorrelation time 

CO 

r^nt = ^ E pW- (69) 

^ t=-oo 

The integrated autocorrelation time controls the statistical error of the sample mean 

/i = TrEO*. (70) 

^^s t=l 
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Indeed its variance can be rewritten for sufficiently large samples as 



-I Ns 1 oo 1 

varif,) = ((A - = j-, C{r - s) c:. - ^ C{t) = — 2r,„,C(0). (71) 

s r,s=l * t=— oo * 

Therefore the 2rj„( is the number of iterations needed to produce one uncorrected 
estimate for O. 

Obviously in real situation the estimation of Tmt have to be done from a finite 
data set. However summation over all available data is misleading since the tail 
of p contains a lot of noise but little signal. Therefore in practice we sum the 
autocorrelations only up to a certain distance Tcut 

Tcut 

TUTcut) = - + Y.Pit)- (72) 
^ t=i 

There are two limitations. If Tcut is too large our estimator is very noisy since its 
variance is 

variriTcut)) = ^-^^^f^r^, (73) 

On the other hand two small value of Tcut introduces a significant bias. The optimal 
value of Tcut can be easily estimated from sufficiently long runs by the "window" 
procedure The recipe is simple: choose Tcut to be the smallest integer such 

that Tcut > CTint{Tcut)- Assuming pure exponential behaviour of p{t) = e"*/^*"' it is 
sufficient to take c = 4. Then the bias of our estimator is of order of = 2% . 

We have calculated the autocorrelation times for two observables: the average 
density of electrons with spin up (ni|) , and that of pairs of electrons with the opposite 
spins on the same site (ni^riii) . The simple updating scheme described above gives 
large autocorrelation times already for K = 1 and f3 = 1. It is important to learn 
whether the long autocorrelations appear while generating fields or during update 
of the A field. We have introduced some modifications to the algorithm to address 
this question. 

One could think that these correlations are caused mainly by the critical slowing 
down introduced by (p fields, especially those with small Pk- This would suggest that 
one would get a better algorithm efficiency by improving updating of the gaussian 
sector. There are many efficient techniques used in simulations of the gaussian fields. 
Among them the Multigrid algorithm eliminates completely critical slowing down 
in the pure gaussian model and therefore we decided to use it to update fields in 



our problem [^. In this approach one considers the sequence of coarse-grid 
problems which approximate original problem on different length scales and the 
local updates of the heat-bath algorithm are supplemented by collective updates. 
We define the set of blocks Bi. For cubic lattice the blocks are taken successively to 
be single sites {Bq ), cubes of side 2 {Bi ), cubes of side 4 {B2) and so on. Obviously 
if the lattice size is not a power of 2 we can choose also other small blocking factor. 
With every block we associate the conditional probability distribution, which for 
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Table 2: 





/3 = 1, f/ = 1 


p = i,u = i 


/3 = 1, f/ = 2 


P=1,U = 2 




lattice 5^ 


lattice 6^8 


lattice 6^8 


lattice 6^14 


Exact determinant 


0.460(2) 


0.473(2) 


0.462(4) 


0.468(5) 




0.2197(2) 


0.2203(4) 


0.195(1) 


0.193(1) 




ri = 3 


ri = 3 


ri = 3 


ri = 2 




T2 = 1.5 


T2 = 1.5 


T2 = 1.5 


r2 = l 


Simple program 


0.461(5) 


0.471(6) 





— 




0.2180(8) 


0.221(1) 








ri = 660 


Ti = 870 








T2 = 320 


T2 = 300 






Simple program 


— 


0.470(4) 





— 


with MG (W-cycle) 




0.2194(6) 










Ti = 160 










T2 = 80 






Global generation 


0.449(7) 








of gaussian fields 


0.221(1) 










n = 180 










T2 = 90 








Simple program 




0.470(5) 


0.463(7) 


0.475(7) 


with preconditioning 




0.220(1) 


0.195(1) 


0.186(3) 






Ti = 100 


ri ~ 200 


Ti ~ 600 






r2 = 60 


T2 ~ 200 


T2 ~ 600 



The results for density of electrons (first hne) and density of pairs {rii^nii) (second 
hne). The autocorrelation times are given for both quantities. Autocorrelation times for 
exact numerical evaluation of the determinant are of the order of unity. 



the 



field reads 



Psif) oc exp 



+ t{xB)^) {Q^Q-akf + f^l 



(74) 



Xb denotes the vector whose components are 1 for sites belonging to the block B 
and elsewhere. The Pb depends only on one real number t. 

This general set up can be directly used in our problem as follows. We begin 
with elementary sites for which we perform the full heat-bath sweep through the 
lattice sequentially updating ^[r] for all i, k and r. Then we organize the lattice into 

(r) 

blocks Bi. We change the 0^ [ at all sites of given block by the same t according to 
distribution ([7^). After having done the sweep over blocks -Bi, for all k and r we 
repeat the same for blocks B2 and so on. This updating scheme is called V-cycle. 
The sweeps can be also organized in such fashion that there are 7' updates of blocks 
Bi. The control parameter 7 defines different classes of MG algorithm. The most 
common choices are 7 = 1 (V-cycle) and 7 = 2 (W-cycle) 

Indeed the MG generation of cf) fields reduces the autocorrelation times substan- 



tially (third row of Table 2). However it introduces additional computational cost. 
In fact the CPU time required to get one independent sample is even bigger than for 
a simple heat-bath method. The maximal decorrelation of auxiliary fields is achieved 
by the independent generation of the eigenmodes. This requires the solution of set 
of linear equations for each field and the results are comparable to the MGMC 
case with W-cycle (4-th row of Table 2). 

Neither simple version nor MG refinement is capable to reproduce the exact de- 
terminant results for U = 2. The system did not thermalize even after 20 times 
more thermalization steps than required for U = 1. This fact can be simply un- 
derstood. When U increases the condition number of M'^M becomes bigger and 
larger number of fields in polynomial approximation is needed. Additionally they 
are stronger coupled to the A field because of the presence of factor exp ^Jup/Nf. 
The constraints imposed by large number of fields on one A field becomes more 
restrictive and the mobility of algorithm rapidly decreases. Performing updates of 
A field more frequently is only a partial solution. 

Better results gives the preconditioning of fermionic matrix. We define the pre- 
conditioning procedure as follows. Let D he a diagonal part of M = A + D. Then 
the matrix M^D~^M is better conditioned than M^M as can be seen from Fig. 3 . 
The effect is clearly visible in the last row of Table 2. As a consequence wider range 
of couplings become available. However the reliable runs at U greater than 2 are 
not feasible. 

The program based on exact evaluation of determinant has no such restrictions 
and works equally good at t/ = 1 like for U = 4. Required CPU time needed 
for both algorithms varies with lattice sizes. On small lattices algorithm based on 
exact evaluation of determinant is substantially faster and it took 20 minutes on 
HP735/125 workstation to obtain the results on lattice 6^8 while comparing with 
10 hours for Liisher implementation. Of course for larger lattices the computational 
effort grows much more faster for exact algorithm. In the region of weak coupling 
and high temperatures we were able to perform simulations on lattices 16^ x 8 with 
the help of the Liischer algorithm. 

The polynomial approximation can be extended simply to the complex plane 
P7| i. e. for nonhermitian matrices. Thus one tries to approximate directly 
the inverse of matrix M. This would reduce the condition numbers for matrices 
entering the problem and would result in the simpler final action. This modification 
of original Liischer idea gave very promising results [0 However in contrast 



to the QCD case the matrix M for the Hubbard Model has eigenvalues with the 
positive as well with negative real parts as can for example be seen on Fig. 4 . 
Because one cannot extend the domain of the applicability of Eq. (Q) beyond the 
singular point (0, 0), it is unfortunately impossible to adapt this modification to the 
Hubbard Model. 

On the other hand we have found that one change the positions of eigenvalues 
by introducing the chemical potential. In fact we managed to shift all eigenvalues of 
the matrix M in simulations of the attractive Hubbard Model to the right half of the 
complex plane. However it occurred at large fi which corresponds to the physically 
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uninteresting filling ((n| + rii) > 0.9). 

Other improvements of the Liischer algorithm have been proposed. Global 



Metropolis step makes the algorithm exact and heat-bath generation of (p fields 
with overrelaxation slightly reduce autocorrelation times (for review see [0). We 
believe however that it would not change the situation qualitatively. 



5.3 Magnetic properties of the Hubbard Model 

In section 2 we argued for the existence of the antiferromagnetic order in the ground 
state of the Hubbard Model. These arguments were based on the mean-field ap- 
proximation and should be independently confirmed in the numerical simulations. 
In fact , the mean field theory is known to overestimate the magnetic ordering and to 
underestimate quantum fluctuations. The situation is rather complicated since the 
Hubbard Model has continuous (spin - SU(2)) symmetry. In that case, in 2c?, the 
Mermin- Wagner theorem prevents the existence of the long-range magnetic 



order at finite temperature. The correlation length can become really infinite only 
at zero temperature when thermal fluctuations vanish. 

Therefore to address the problem of magnetic properties of the Hubbard Model 
it is necessary to perform simulations at large (3. This in turn requires sufficient 
Nt in the discretized form of the path integral formulation to reduce the systematic 
errors of order of {P/Nt)"^. Moreover in our particular representation the number 
of electrons with spin down and with spin up differ. It is special inconvenience in 
measurements some correlation functions leading to the errors which are extensive 
function of the lattice size. 

To estimate the temperature we can reach, in the physical units we assume that 
8K = leV. Then it follows from Eq. (|^) that j3 = 5 and K = 1, the most typical 
values for our runs correspond to temperature T = -^eV = 290 Kelvin degrees. It 
would be difficult to reach much lower temperatures due to computer restrictions 
and finite Nt effects. 

We begin discussion with a local quantity which can be simply measured on 
small lattices with good accuracy. The average number of pairs with opposite spins 
(n^rii) on the same site shows how repulsive interaction forbids double occupancy of 
a site. Clearly increasing U leads to reduction in number of pairs. Fig. 5 compares 
results from QMC simulations with the mean field prediction. The derivation of the 
latter is presented in the Appendix A. Here we note a good qualitative agreement 
with MC results. For this quantity a 4 x 4 lattice gives rather good approximation 
to the bulk result and we did not notice meaningful finite size effects as can be seen 
from Table. 3 . 

To decide weather the ground state has the long range antiferromagnetic order 
one can attempt to measure the equal time correlations between magnetization on 
different sites 

C{lx, ly) = {{rii^ - nii){ni+^ - Ui+n)). (75) 
This quantity obtained form simulations on 4 x 4 lattice at j3 = 5 and U = 4is shown 
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{nirriii) 


N=16 


N=36 


N=64 


N=100 


N=144 


U = l(3 = lNt = 8 


0.2195(4) 


0.2204(4) 


0.2202(5) 


0.2202(3) 


0.2198(3) 


U = 4f3 = 5Nt = 30 


0.134(2) 


0.138(2) 









Table 3: The density of pairs on the same site as a function of the spatial volume of the 
lattice = N^. For this quantity the lattice 4^ already gives good approximation to the 
bulk result 



on Fig. 6 . The visible zig-zag is an indication of the onset of the antiferromagnetism. 
We found that this observable is very sensitive to finite Nt effects. We needed 60 
time slices to obtain clear signal (the analogical zig-zags for A^^^ = 30 and Nt = 40 
do not look so nice). This is not surprising since the inequality between (rii^) and 
(nj|) is equivalent to putting the system into a magnetic field which destroys the 
antiferromagnetic correlations. Within the spin wave theory the correlation between 
two most distant points on the lattice is simply related to the antiferromagnetic order 



parameter m [45 



777^ 

C{NJ2, = — + 0(l/iV,) (76) 

with the finite size corrections of order 0{1/Nx) on the lattice with spatial volume 
= N"^. For the parameters of run the mean field result m^/3 = 0.15 should be 
compared with the value 0.10 ±0.04 read from Fig. 6 . Thus even for relatively high 
temperature of our simulations we obtained a good agreement. We did not estimate 
the finite size effects because it requires the simulations on bigger lattices. 

Fourier transformation of C{1) gives the magnetic structure factor 

S{q) = Y: exp {tql)C{l) (77) 

which is especially well suited for detecting excitations with wave vector q. In the 
presence of antiferromagnetic long-range order we expect divergence of the 5'(7r, vr) 
with the lattice size according to 

2 

7T? 

S{n,n)=N— + Oil/N,). (78) 

So far no one has shown in direct simulations the existence of antiferromagnetic 
phase at U below 2 |Q. Since we can reach quite large lattices at small U using the 
Liischer method we made such an attempt. Fig. 7 gives the S{kx, k^) as a function 
of momentum on lattices with spatial sizes ranging from 4 x 4 to 16 x 16 for /3 = 1 
and U = 1. Within the statistical errors there is no dependence of the magnetic 
structure factor on lattice volume at this temperature. An exception is S'(0, 0) which 
rises with almost linearly. It is a simple artefact of finite Nt and does not mean 
that system has ferromagnetic long range order. Indeed the S'(0, 0) is defined as 
an extensive quantity and naively we can expect the contribution to the S'(0, 0) of 
order of A^(An)^ coming from the difference An between and (n^). This crude 
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Nt = 6 


Nt = 8 


iVi = 12 


a = 


0.33(3) 


0.34(2) 


0.34(3) 


b= 


0.044(4) 


0.030(2) 


0.013(2) 


{Anf 


0.041 


0.025 


0.010 



Table 4: Parameters a and b obtained from simulations on lattices with Nt = 6, 8, 12. 
The mean difference (An)^ is also given. 



N=16 


N=36 


2.6 ±0.2 


3.2 ± 1.0 


2.8 


4.2 



Table 5: Comparison of results for S'(7r,7r) obtained from our studies (first row) with 
those from Ref. [12] (second row). 

estimation is tested in Fig. 8 where the S{0, 0) is shown as a function of the spatial 
volume of the lattice for Nt = 6, Nt = 8 and A^^ = 12. For each A^^ we found the 
parameters of linear fit 

S{0,0) = bN + a (79) 

The results are listed in Table 4. The values of b should be compared with the 
mean difference (An)^. The remaining contribution a, which does not depend on 
Nt may be regarded as a true 5(0, 0). The nonozero value of S{0, 0) is an effect of 
the short-range correlations. 

The simulations described above were performed at rather high temperature 
(corresponding j3 = 1). We did not found the antiferromagnetism because of the 
disordering effect of thermal fluctuations. We performed some simulations with the 
help of the exact algorithm at [/ = 4 and the /3 = 5 where the scaling predicted 
by the Eq. (^) has been observed [|12]. Due to strong finite Nt effects to obtain 
valid result one has to perform extrapolation. The finite Nt errors introduced by 
our lattice discretization are of the order of {-^Y- In Fig. 9 we present ^(Tr, vr) as a 
function of parameter {^Y- In the limit 1/Nt for the lattice 4 x 4 we obtained 
2.6 ± 0.2 and for the lattice 6x6 3.2 ± 1.0. They compare satisfactorily with MC 
simulations of Ref. []12[ (see Table. 5). 

5.4 Fermi surface 

It is also very interesting to study information contained in one particle finite tem- 
perature Green's function defined in the momentum space and imaginary time as 

G{k,T) = -{Tak.{T)aU0)), (80) 

where a(r) is the annihilation operator in the Heisenberg picture and T time or- 
dering operator. The limit of G{k, r 0~) has simple physical content. It gives 
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distribution of electrons in the momentum space. To see the shape of the Fermi sur- 
face one has to perform simulations on sufficiently large lattices. The largest lattice 
accessible to our simulations was 8^ x 20. For the parameters of run U = 4 and 
j3 = 5 the number of time slices might appear too small. However, for this quantity 
we do not expect the presence of the mechanism described above which leads to the 
errors which are extensive function of the spatial lattice volume. 

Finally, one can also compare our results with the mean-field prediction. It 
follows from Eq. (^) in Appendix A that 




The energy gap A = 1.38 is almost the same for the lattice 6x6 and for the lattice 
8x8. Consequently we can together compare results from lattices 6x6 and 8x8 
with the mean field formula (solid line in Fig. 10 ). Even for our relatively large 
P/Nt = 0.25 the results we obtained are in a good agreement with the theoretical 
prediction. However the Fermi surface of the interacting system is more smooth 
than in the noninteracting case (dotted line). Obviously this effect is an increasing 
function of U. 

The lack of sharp Fermi surface suggests vanishing Z and behaviour different 
from the Fermi liquid. Of course it would be very interesting to measure the n{k) at 
non zero chemical potential fi and find the critical value Hcrit when the sharp Fermi 
surface appears. Unfortunately the QMC simulations are not conclusive at small /i 
due to sign problem. Currently it is only confirmed that ficrit is small. 



5.5 Effective hopping 

As pointed out in various approximate schemes the interaction leads to the reduction 
in the effective hopping. In our Monte Carlo simulations we measured the ratio 



'U 



(82) 



which shows how interaction reduces the effective hopping matrix elements between 
nearest neighbor sites This quantity can be measured directly or can be ex- 
pressed in terms of Green's function. Indeed Keff can be written as an expectation 
value over the whole lattice 

Keff = j^{Yl O'Ux+u + h.c), (83) 

x,u 

where x denotes the sites on the lattice, u is the unit versor and N is spatial volume 
of the lattice. For simplicity we drop the spin index. Replacing the operators by 
their Fourier transforms yields 

Keff = E(4«.)(e^'^'^ + e-^'=''^)Ee^^'"''^^ 

k,k'v a; 



M) 
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Finally 

Keff = -^Y.<k)e,. (85) 

k 

The last equation is very convenient since the mean field formula for n{k) has been 
already given. In Fig. 11 we compare the results from simulations and the mean-field 
prediction. The observed reduction is smaller than given by the theory. It is not 
surprising since the mean field approximation is known to overestimate effect of U . 
Indeed the mean field ground state does not take into account the suppression of 
double occupancy. Consequently the effect of the interaction term is overestimated. 



6 Summary and conclusions 

Motivated by the interest for numerical approach proposed by Liischer we have de- 
veloped a Monte Carlo algorithm for studying the two dimensional Hubbard Model. 
Our method is based on the field theoretical formulation in (2 + 1) dimensions, to 
which we applied the polynomial approximation proposed by Liischer 0. Finally 
we obtained the local action and therefore we were able to study the system using 
bosonic techniques. We have introduced a number of improvements which reduce 
the long autocorrelation times. In particular we show that by introducing a pre- 
conditioning of the fermionic speeds up the algorithm significantly. We have also 
compared the efficiency of the Liischer algorithm with the simple algorithm based 
on the direct update of the determinant of the fermionic matrix. It follows from this 
comparison that at present stage Liischer technique does not seem to be a true alter- 
native to the standard simulations of the Hubbard Model. The Liischer method can 
not reach the most interesting region of strong coupling because then the fermionic 
matrices are very badly conditioned, much more worse than in the case of QCD. A 
large number of the additional bosonic fields is required and the algorithm suffers 
from strong autocorrelations. Another difficulty arises from nonhermicity of the 
fermionic matrix of the Hubbard Model which complicates the resulting action and 
squares the condition numbers for matrices entering the problem. The nonhermitian 
version of the Liischer algorithm can not be also applied. 

The algorithm based on the direct update of the determinant of the fermionic 



matrix was modeled after the algortihm from Ref. [|TT|. However we found that nu- 
merical instabilities do not appear in our particular representation of path integrals, 
which allow to avoid additional complications. Hence we managed to reach region 
of strong coupling and quite large {3 at half-filling with the help of this very simple 
algorithm. Whenever it was possible we also compared our results with results of 
previous investigations as well with the mean field theory. 

To show the that the ground state has antiferromagnetic long range order we 
measured the correlation between two different sites on the lattice. We found that 
main limitations in measurements of this quantity are finite Nt effects which intro- 
duce significant bias. Nevertheless, through the extrapolation iVt — >■ oo, we were 
able to obtain results which are in a qualitative agreement with other QMC simu- 



lations []TT| |T2[. We have also attempted to measure the one-particle properties. In 
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particular we showed that the Fermi surface is not sharp with the shape described 
well by the mean-field theory. This shows that the system is in the insulating Neel 
state. Thus we did not found any unexpected behavior and our results confirm the 
often presented opinion that the physics of the half-filled Hubbard Model is well 
understood. 

On the other hand there is the case of non half-filled band with possible exis- 
tence of the superconducting phase. Although the Projector Quantum Monte Carlo 
simulations delivered substantial progress in this field the situation still remains 
unclear. Especially just below the half-filling where the sign problem is a serious 
obstacle. Current simulations are consistent with the most orthodox point of view 
treating holes as quasiparticles but other non Fermi liquid scenarios are not com- 
pletely excluded fl^. It is also impossible to make final conclusions about 
existence of the superconducting phase because the range of parameters accessible 
to simulations is too narrow. 
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A Appendix A 



In this Appendix we apply the mean field theory to the half-filled Hubbard 



Model on the square lattice with the antiferromagnetic order in the ground state. 
Within this approximation one rewrites the operator of the number of electrons at 
a given site i as riia = {uia) + {niu — {ni„)) where {ni„) is the expectation value in 
the ground state. Then assuming the fiuctuations to be small, the interaction term 
in Hamiltonian becomes 

ni^Uii = -{ni^){nii) + rii^iriii) + ^^^(ni-f) (86) 

It is convenient to choose as the ground state Spin Density Wave (SDW) state 
defined by 



1 



1 r 



where as usual (—1)* = (—1)*^+*!' is the parity of site, and m is a variational param- 
eter. Its value will be determined from the minimum energy condition. Note that 
we are able to reproduce the result for interacting electrons putting m = and the 
antiferromagnetic Neel state (m = 1). Transformation of the original Hamiltonian 
to the momentum space yields the mean field Hamiltonian in the form 

fee k 

where Cfc = — 2_ft'(cos + cos ky) are the bare energy levels and Q = (tt, vr) . The 
above Hamiltonian is diagonal in the spin indices and the operators with momentum 
k can interact only with those of momentum k + Q. Thus our problem reduces to 
the diagonalization of the set of matrices 



^9) 



/ efc ±mU/2 
\ ±mU/2 tk+Q 

The eigenvectors and eigenvalues can now be simply obtained 

Tfct ^ = ^fcOfcT + Ukttk+Qh iti'^ = ^kaki - Ukak+Qi- (90) 

There are two eigenvectors per spin and momentum k is in the reduced Brillouin 
zone. The transformation amplitudes Uk and Vk are 

2 



k 



A = -^. (91) 



where A is the energy gap parameter and ±Ek are the eigenvalues of 7^^''. At the 
edge of the Fermi surface k = {— 71/2,71/2) we have two possible energies E = ±A 
thus the system acquires gap equal 2 A. The mean field ground state is obtained by 
filling the states with lower energy (i.e. those with the sign (-) in front of Ek). 

\^mf)= n (7&Y(7ilV|0), (92) 

\k\<kp 

with the empty state denoted by |0). The condition for energy minimum leads to 
the self-consistent equation for the energy gap 

U ^ N ^ ^^^^ 
\k\<\kp\ (e/c + ^ J 

where the sum over momenta is restricted to half of Brillouin zone {kp is the mo- 
mentum at the noninteracting Fermi surface). The above equation for A for finite 
lattices can be easily solved numerically. The mean field predictions can be simply 
obtained for physically interesting quantities. For example the double occupancy of 
the site is simply related to the antiferromagnetic order parameter m 

{rii^rii^) = {ni^){n,i) = ^(1 - m^). (94) 
One can also obtain the occupation number in the momentum space 



rika = {■ipMFWkaO'kal'ipMF) = 2 V''" " 'E~ ) ' ^^^^ 

Indeed the electrons with momentum k are created from the state |0) with the 
amplitude Vk, given by Eq. (0) if the momentum k lies inside the Fermi surface. 
Similarly the electrons outside the Fermi surface are created with the amplitude 

Uk-Q = Vk- 



B Appendix B 



In this Appendix we review the formalism of Grassmann variables [|T7| |T^]. A set 
of anticommuting Grassmann variables is defined by the relation 

[ipi, ipj]^ = + = 0. (96) 

For each ipi there is also corresponding independent variable tp*. Moreover we have 

{^l-'-^nY = K-'-^l ■ (97) 

Assuming the properties of linearity and invariance under translation of variables 
the integration can be defined through the formulas 



dipip = i 
J dipl = 

dip*dtptp*ip = 1. (9J 



Since the general function of Grassmann variables is in fact polynomial the above 
equations are sufficient to compute any integral. In addition to integration we 
consider the differentiation with the respect to the Grassmann variable. This can 
be defined by the action on a constant function and an anticommutation relation 



= 1. (99) 

+ 



Function of Grassmann variables form the Hilbert space on which transfer matrix 
acts. To make this statement more specific we consider as an example the case of 
one degree of freedom. Then one can relate the function f{ip*) = 1 to the vacuum 
state |0) and the function /(t/'*) = V'* to the occupied state The operators 
and a would correspond to 'ijj* and ^ respectively. 

It is convenient to introduce the simple analog of Dirac S function for the anti- 
commuting variables 

5(^*,^'^) = / (rf^)eE(V'*-V''*)V'^ (100) 



where (dip) denotes Yidipi in some prescribed order. Note that the 5 function has 
the expected property 

fir)^ Js{r,i^'*WTf{i^'*). (101) 

Now it is easy to verify the correspondence on the level of integral representations 
for operators acting on a general state |/) 

^ / (ci^)eS(^*-^'*)^ (#T/(V''*) ■ (102) 

The result is very simple and states that under integral a and should be replaced 
by ip and ip* respectively. This immediately generalizes to any normal ordered 
function of a and 

: g(a\ a) : I/) - / gi^, V') (#)eE(^^-^'*)^ (#T/(V''*)- (103) 

It is very important that the function g is normal ordered. Variables ip and ip* simply 
anticommute while the anticommutation relation between tjj* and introduces 
additional factors. 



The trace of a normal ordered operator can be rewritten in the form of integral 
over anticommuting variables 

Tr : g{a\a) : = J {d^yg{^\^)e^^^*^{d^). (104) 

A little algebra yields the trace of normal ordered factors as a multiple integral over 
anticommuting variables 

Tr[: gi{a\ a) :: g2{a\ a) : ■ ■ ■ : gNt{a\ a) :] = 

/ n[(#t)*^t(^^^t)(#t)eS'^?(^-^*-)], (105) 



with antiperiodic boundary conditions ipo = —ipNf For operators g being the trans- 
fer matrices of the Hubbard Model the integration is performed by the standard 
formula 

J {dipf {d^lj)e^*^^^ = det M. (106) 

which directly gives Eq. (|23|) . 

The expectation value of observables just inserts another factor in the above 
integral. The additional factors under integral are 

j {d^Nt+lT (d^Nt+l) /i(^^,+l, ^iVi+l) exp [tpN.+li^Nt+l - ^Nt) + i'lii^Nt+l - ^TvJ]- 

(107) 

As an example we consider the expectation value for the numbers of electrons 
with spin up and the corresponding function : h[a\ a) := a^a. Then the integration 
over additional factors reduces simply to the insertion 1 — iplipNt under the main 
integral. The evaluation of integral of such type is simple if we make use of the 
identity 

j (#)*#e^**^^+^^+'^*'' = det M e"^^"'^. (108) 

The additional factors can be obtained simply by differentiation with the respect to 
the anticommuting variables t] and ^ next putting them zero. In our example 

J (#)e^*^^ (1 - i^l,^^,N.) = det M (1 - M-^.^^^,). (109) 

Cyclically shifting in t direction yields 

(n,t) = (l + Mr;|.,^i). (110) 

In similar derivation for the electrons with spin spin down one should first perform 
the particle hole transformation aja^ = bb^ = 1 -- Wb. The corresponding insertion 
is i'l'ipNt and the final result reads 

(n.i) = -(Mr,i,+,). (Ill) 

Similar considerations give the expressions for the other observables of interest 

(^.T^^i) = -(M^i,m(l + ^i,m)) (112) 

C(z,j) = (1 + 2Mr,i,+, + 2Mr,;^.,+, - 2M,r;i A. 

+AM-^,^,_,,M7^^^,^, - 2M^|^.,+iMr,;,,+i) (113) 

M^) = lT.e''^''-'HM-i.,y,^,-{-irhirM (114) 

xy 
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Figure Captions 



1. Fermi surface for the 2d Hubbard Model on the square lattice. The dashed 
area corresponds to the occupied states for half-fiUing. 

2. Errors of the Liischer approximation measured on a single typical configuration 
genrated at U — 1 and P — 1. The quantity 5 defined in text is shown as a 
function of e. The solid, dashed and dotted line is for number of fields 50, 70 
and 100 respectively. 

3. The distribution of the condition numbers for the M^M matrix (solid line) and 
the preconditioned matrix M^D^^M (dashed line). Results from simulations 
on 6^ X 12 lattice at [/ = 4 and /3 = 1. 

4. Eigenvalues of the matrix M plotted on the complex plane. A typical config- 
uration was taken from simulations at U — 1 and (3 — 1. 

5. The double occupancy versus U on a 4 x 4 lattice with P — 5. The solid hne 
is the mean-field result. 

6. Magnetic correlation function C{1) on the 4x4 lattice. The point / traces 
the triangular path shown in the picture. The antiferromagnetic correlations 

are clearly visible at Nt = 60, j3 = 5 and U = 4. The shape of curve is very 
sensitive to the finite Nt effects. The results at Nt — 30 (dotted line) and 
Nt = 40 (dashed line) are also shown for comparison. 

7. Magnetic factor S(kx, kx) along the main diagonal of the Brillouin zone. Re- 
sults from lattices 6x6 (empty circles), 8x8 (triangles), 12 x 12 (circles) 
and 16 x 16 (crosses) are collected. Other parameters of runs were U = 1 
and (3 — 1. There is no signal of long range antiferromagnetic order at this 
temperature. 

8. Finite Nt analysis of antiferromagnetic factor S{7r, vr). Its value is extrapolated 
in the parameter (-^) . Results from lattices 4x4 (circles) and 6x6 (triangles) 
are presented. The extrapolation with A^^ — > oo gives result 2.6 ± 0.2, and 
3.2 ± 1.0 respectively. 

9. The ferromagnetic factor S{Q, 0) as a function of spatial volume of the lattice. 
Three sets of data are shown for Nt — 6 (circles), Nt — S (triangles), and 
At = 12 (crosses). 

10. The distribution of electrons {nk^ + Uki) in the momentum space. Monte Carlo 
results from lattice 6x6 (crosses) and 8x8 (triangles) are shown. They are well 
described by the mean-field theory (solid line). Dashed line is a noninteracting 
Fermi distribution /(e^) — 2{exp{—(3ek) + 1)~^- 

11. The effective hopping versus U on a 4 x 4 lattice with (3 — b. The sohd hne is 
the mean field result. 
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